tutorial_scl - Seurat_Markdown

suppressPackageStartupMessages(library(scater))
suppressPackageStartupMessages(library(Seurat))

data_path = "D:/Dropbox/매뉴얼/SC_tutorial"
figure_path = "D:/Dropbox/매뉴얼/SC_tutorial"
setwd(data_path)

Load scesetFiltered, hvg

scesetFiltered = readRDS( file= "scsetFiltered_norm.rds")
load(file="hvg.RData")

Create seurat object

seuset <- CreateSeuratObject(
  raw.data = counts(scesetFiltered)
)
seuset@data = exprs(scesetFiltered)
seuset@scale.data = exprs(scesetFiltered)
seuset@var.genes = rownames(hvg)

PCA

PCA = 50
seuset <- RunPCA(seuset, pcs.compute = PCA, weight.by.var = FALSE, do.print = TRUE)

## [1] "PC1"
##  [1] "Hbb-bt"  "Opalin"  "Vtn"     "Mog"     "Trem2"   "Mag"     "Selplg" 
##  [8] "Aif1"    "Mobp"    "Gsn"     "Mal"     "Fcrls"   "Tmem119" "Fcer1g" 
## [15] "Csf1r"   "Ctla2a"  "Rgs5"    "Flt1"    "Cldn11"  "Hbb-bs"  "Pglyrp1"
## [22] "Hmgb2"   "Ptgds"   "Tyrobp"  "Fgfr3"   "Ly6c1"   "Fam107a" "Cldn10" 
## [29] "P2ry12"  "Id1"    
## [1] ""
##  [1] "Malat1"  "Fth1"    "Tmsb4x"  "mt-Nd1"  "Calm2"   "mt-Cytb" "Ftl1"   
##  [8] "Sepw1"   "Snca"    "Rps5"    "Rps9"    "Rplp0"   "Ubb"     "Pcp4"   
## [15] "Rps4x"   "Ppp3ca"  "Meg3"    "Rpl10"   "Ptma"    "Aldoa"   "Ybx1"   
## [22] "Tuba1a"  "Zbtb20"  "Bex2"    "Tmsb10"  "mt-Nd2"  "Mllt11"  "Ypel3"  
## [29] "Mt3"     "Rps29"  
## [1] ""
## [1] ""
## [1] "PC2"
##  [1] "Cst3"     "Apoe"     "Dbi"      "Cd81"     "Ptn"      "Aldoc"   
##  [7] "Glul"     "Atp1a2"   "Slc1a2"   "Fabp7"    "Mt1"      "Clu"     
## [13] "Cpe"      "Ftl1"     "Itm2b"    "Sparcl1"  "Pantr1"   "Tsc22d4" 
## [19] "Slc1a3"   "Ndrg2"    "Cspg5"    "Ckb"      "Lamp1"    "Scd2"    
## [25] "S100a16"  "Car2"     "Plpp3"    "Prdx6"    "Cd63"     "Serpine2"
## [1] ""
##  [1] "Snca"   "Pcp4"   "Meg3"   "Ppp3ca" "Calm2"  "Ncdn"   "Olfm1" 
##  [8] "Camk2b" "Bex2"   "Sncb"   "Mllt11" "Cplx2"  "Cnih2"  "Snrpn" 
## [15] "Ndrg4"  "Ctxn1"  "Chn1"   "Map1b"  "Snap25" "Ly6h"   "Ywhah" 
## [22] "Camk2a" "Ncald"  "Plppr4" "Rtn1"   "Atp1b1" "Grin2b" "Nsf"   
## [29] "Caly"   "Atp2b1"
## [1] ""
## [1] ""
## [1] "PC3"
##  [1] "Mt3"     "Mt1"     "Slc1a2"  "Tspan7"  "Chchd10" "Ncdn"    "Atp1b1" 
##  [8] "Sparcl1" "Aldoc"   "Ppp3ca"  "Aldoa"   "Chn1"    "Apoe"    "Ppp1r1a"
## [15] "Gpm6a"   "Clu"     "Camk2a"  "Atp1a2"  "Malat1"  "Fth1"    "Camk2b" 
## [22] "Cplx2"   "Adcy1"   "Meg3"    "Ndrg4"   "Camk2n1" "Ttyh1"   "Olfm1"  
## [29] "Nsf"     "Plppr4" 
## [1] ""
##  [1] "Tmsb10"   "Tuba1a"   "Tubb5"    "Ftl1"     "Rplp0"    "Ptma"    
##  [7] "Igfbpl1"  "Rps5"     "Tubb2b"   "Marcksl1" "Rps9"     "Tmsb4x"  
## [13] "Ybx1"     "Rps4x"    "Eef2"     "Tubb3"    "Rpl10"    "Nnat"    
## [19] "Stmn2"    "Hn1"      "Ubb"      "Sox11"    "Prdx2"    "Junb"    
## [25] "Rps28"    "Sox4"     "Calb2"    "Xist"     "Gm1673"   "Fxyd6"   
## [1] ""
## [1] ""
## [1] "PC4"
##  [1] "Tubb2b"        "Tuba1a"        "Zbtb20"        "Tmsb10"       
##  [5] "Nnat"          "Igfbpl1"       "Tubb3"         "Stmn2"        
##  [9] "Xist"          "Marcksl1"      "Tubb5"         "Rtn1"         
## [13] "Fabp5"         "Gm1673"        "Slc1a2"        "Mllt11"       
## [17] "Hn1"           "Mt3"           "Cnih2"         "Sox11"        
## [21] "Aldoc"         "Cpe"           "6330403K07Rik" "Fxyd6"        
## [25] "Prdx2"         "Mmd2"          "Calb2"         "Pantr1"       
## [29] "Clu"           "Cspg5"        
## [1] ""
##  [1] "Sparc"  "Bsg"    "Sepp1"  "Tmsb4x" "Junb"   "Arpc1b" "Hexb"  
##  [8] "Itm2b"  "C1qa"   "Fth1"   "Lgmn"   "C1qb"   "C1qc"   "Ctss"  
## [15] "Ctsd"   "Cldn5"  "Cst3"   "Jun"    "Gatm"   "Cd81"   "P2ry12"
## [22] "Cyba"   "Malat1" "B2m"    "Tyrobp" "Klf2"   "Trf"    "Itm2a" 
## [29] "Igfbp7" "Rgs10" 
## [1] ""
## [1] ""
## [1] "PC5"
##  [1] "Cst3"    "C1qa"    "Hexb"    "C1qb"    "C1qc"    "Ctss"    "Ctsd"   
##  [8] "P2ry12"  "Lgmn"    "Tyrobp"  "Rgs10"   "Fcrls"   "Tmem119" "Csf1r"  
## [15] "Fcer1g"  "Selplg"  "Aif1"    "Junb"    "Apoe"    "Trem2"   "Sirpa"  
## [22] "Rps4x"   "Ctsz"    "Rps29"   "Mt3"     "Ctsb"    "Rpl10"   "Slc1a2" 
## [29] "Rplp0"   "Rps28"  
## [1] ""
##  [1] "Plp1"    "Bsg"     "Mbp"     "Cnp"     "Sept4"   "Cldn5"   "Cryab"  
##  [8] "Car2"    "Gng11"   "Tsc22d1" "Ptn"     "Cldn11"  "Mal"     "Itm2a"  
## [15] "Mobp"    "Igfbp7"  "Ptgds"   "S100a16" "Mog"     "Qdpr"    "Aplp1"  
## [22] "Pllp"    "Malat1"  "Mag"     "Ly6c1"   "Ptma"    "Opalin"  "Fth1"   
## [29] "Olig1"   "Car4"   
## [1] ""
## [1] ""

qplot(x = seq(1:PCA),y = seuset@dr$pca@sdev,
      xlab = "PC", ylab = "Eigenvalue")

PC25까지 사용하여 TSNE 및 Clustering

PCA_use = 25
seuset <- RunTSNE(seuset, dims.use = 1:PCA_use, do.fast = T, seed.use = 123456, perplexity=100)
seuset <- FindClusters(seuset, reduction.type="pca", dims.use = 1:PCA_use, save.SNN = TRUE, force.recalc = TRUE)

# save(seuset, file= "seuset.RData")
load("seuset.RData")

clusters_seurat = seuset@ident
# save(clusters_seurat, file= "clusters_seurat.RData")
load("clusters_seurat.RData")

Seurat 패키지 function으로 plot하기

PCAPlot(seuset)

TSNEPlot(seuset, do.label = TRUE)

ggplot2 사용하여 plot하기

# publication celltype labels
qplot(x=seuset@dr$tsne@cell.embeddings[, "tSNE_1"],
      y=seuset@dr$tsne@cell.embeddings[, "tSNE_2"],
      xlab ="component1", ylab = "component2",
      colour = scesetFiltered$celltype)

각 cluster별 Markers 찾기

# markers <- FindAllMarkers(
#   object = seuset, 
#   test.use = "wilcox", 
#   only.pos = TRUE, 
#   min.pct = 0.25, 
#   logfc.threshold = 0.25
# )

# save(markers, file="markers.RData")
load("markers.RData")

상위 5 markers heatmap 그리기

suppressPackageStartupMessages(library(dplyr))
top5 <- subset(markers, cluster %in% c(0,1,2,3,4,5,8,10)) %>% group_by(cluster) %>% do(head(.,5))

DoHeatmap(
  object = seuset,
  cells.use = names(seuset@ident)[seuset@ident %in% c(0,1,2,3,4,5,8,10)],
  genes.use = top5$gene, 
  slim.col.label = TRUE, 
  remove.key = TRUE
)

Seurat 기능

RidgePlot(seuset, features.plot = c("Calb1"))

## Picking joint bandwidth of 0.0474

features.plot = c("Cdk1", "Ascl1", "Tfap2c", "Eomes",
                  "Igfbpl1","Calb2","Plk5",
                  "Gfap","Hes5", "Sox2", "Emx2",
                  "Foxg1", "Egfr", "Prom1",
                  "Lpar1","Nes", "Neurog2","Top2a",
                  "Mcm2", "Dcx", "Neurod1", "Calb1"
                  )
DotPlot(object = seuset, genes.plot = features.plot, plot.legend = TRUE, x.lab.rot =T)

DotPlot(object = seuset, genes.plot = unique(top5$gene), plot.legend = TRUE, x.lab.rot =T)

FeaturePlot(object = seuset, features.plot = "Igfbpl1", no.legend = FALSE,
            do.hover = TRUE, data.hover = c("ident", "PC1", "nGene")
            , dark.theme = TRUE)